﻿#include "calcutils.h"
#include <QtGlobal>
#include <QStringList>

static inline bool isFinite(double v)
{
    return std::isfinite(v);
}

bool CalcUtils::toScaled(double v, AxisScale s, double &out, QString *err)
{
    if (s == AxisScale::Linear) {
        out = v;
        return true;
    }
    if (v > 0.0) {
        out = std::log10(v);
        return true;
    }
    if (err)
        *err = "Value <= 0 used on logarithmic axis.";
    out = 0.0;
    return false;
}

double CalcUtils::fromScaled(double v, AxisScale s)
{
    return (s == AxisScale::Linear) ? v : std::pow(10.0, v);
}

// ====== 拟合/插值 ======
FitResult CalcUtils::fitLine(const QPointF &p1, const QPointF &p2,
                             AxisScale xScale, AxisScale yScale)
{
    FitResult r;
    double x1p = 0, y1p = 0, x2p = 0, y2p = 0;
    QString msg;
    if (!toScaled(p1.x(), xScale, x1p, &msg)) {
        r.message = "p1.x: " + msg;
        return r;
    }
    if (!toScaled(p2.x(), xScale, x2p, &msg)) {
        r.message = "p2.x: " + msg;
        return r;
    }
    if (!toScaled(p1.y(), yScale, y1p, &msg)) {
        r.message = "p1.y: " + msg;
        return r;
    }
    if (!toScaled(p2.y(), yScale, y2p, &msg)) {
        r.message = "p2.y: " + msg;
        return r;
    }

    const double dx = x2p - x1p;
    if (dx == 0.0) {
        r.message = "x1' == x2' (vertical in scaled domain)";
        return r;
    }

    r.k = (y2p - y1p) / dx;
    r.b = y1p - r.k * x1p;
    r.ok = true;
    return r;
}

FitResult CalcUtils::interpolateGivenX(const QPointF &p1, const QPointF &p2,
                                       double x3, double &y3Out,
                                       AxisScale xScale, AxisScale yScale)
{
    FitResult r = fitLine(p1, p2, xScale, yScale);
    if (!r.ok)
        return r;
    double x3p = 0;
    if (!toScaled(x3, xScale, x3p, &r.message)) {
        r.ok = false;
        return r;
    }
    const double y3p = r.k * x3p + r.b;
    y3Out = fromScaled(y3p, yScale);
    return r;
}

FitResult CalcUtils::interpolateGivenY(const QPointF &p1, const QPointF &p2,
                                       double y3, double &x3Out,
                                       AxisScale xScale, AxisScale yScale)
{
    FitResult r = fitLine(p1, p2, xScale, yScale);
    if (!r.ok)
        return r;
    double y3p = 0;
    if (!toScaled(y3, yScale, y3p, &r.message)) {
        r.ok = false;
        return r;
    }
    if (r.k == 0.0) {
        // 水平线：若 y3 不在直线上则无解；等于则任意 x，这里返回 p1.x
        const double yConst = r.b;
        if (std::fabs(y3p - yConst) > 1e-12) {
            r.ok = false;
            r.message = "k==0 and y3 not on line";
            return r;
        }
        x3Out = p1.x();
        return r;
    }
    const double x3p = (y3p - r.b) / r.k;
    if (!isFinite(x3p)) {
        r.ok = false;
        r.message = "x3' non-finite";
        return r;
    }
    x3Out = fromScaled(x3p, xScale);
    return r;
}

FitResult CalcUtils::fitAndInterpolate(const QPointF &p1, const QPointF &p2,
                                       AxisScale xScale, AxisScale yScale,
                                       const double *x3Given, const double *y3Given,
                                       double *x3Out, double *y3Out)
{
    FitResult r = fitLine(p1, p2, xScale, yScale);
    if (!r.ok)
        return r;
    const bool hasX = (x3Given != nullptr);
    const bool hasY = (y3Given != nullptr);
    if (hasX && !hasY) {
        double y;
        auto r2 = interpolateGivenX(p1, p2, *x3Given, y, xScale, yScale);
        if (!r2.ok)
            return r2;
        if (y3Out)
            *y3Out = y;
        return r2;
    } else if (!hasX && hasY) {
        double x;
        auto r2 = interpolateGivenY(p1, p2, *y3Given, x, xScale, yScale);
        if (!r2.ok)
            return r2;
        if (x3Out)
            *x3Out = x;
        return r2;
    }
    return r;
}

/**
 * @brief 计算在给定频率范围内扫描所产生的点
 * 该函数精确模拟了逐点生成的逻辑，以确保与实际测试设备的行为完全一致，
 * 特别是对于最后一个点的处理（当计算值超过停止频率时，将停止频率作为最后一个点）。
 * @param freqStart 起始频率。必须 > 0。
 * @param freqStop  终止频率。必须 >= freqStart。
 * @param stepSize  步进值。
 *                  - 在 LINEAR 模式下，这是绝对的频率增量（例如 1000.0 代表 1000Hz）。必须 > 0。
 *                  - 在 LOGARITHMIC 模式下，这是百分比（例如 5.0 代表 5%）。必须 > 0。
 * @param mode      步进模式（StepMode::LINEAR 或 StepMode::LOGARITHMIC）。
 *
 * @return long long 计算出的总点数。
 * @throw std::invalid_argument 如果输入参数无效。
 */
long long CalcUtils::calculateFrequencyPoints(double freqStart, double freqStop, double stepSize, AxisScale mode)
{
    // 参数有效性校验
    if (freqStart <= 0 || freqStop < freqStart || stepSize <= 0) {
        throw std::invalid_argument("Invalid arguments: Ensure freqStart > 0, freqStop >= freqStart, and stepSize > 0.");
    }

    // 如果起点和终点几乎相同，则只有一个点
    // 使用一个小的容差值来安全地比较浮点数
    if (std::abs(freqStart - freqStop) < 1e-9) {
        return 1;
    }

    long long pointCount = 0;
    double currentFreq = freqStart;
    bool stopFreqReached = false;

    // 模拟循环来计数
    while (!stopFreqReached) {
        pointCount++;

        if (std::abs(currentFreq - freqStop) < 1e-9) {
            stopFreqReached = true;
            continue;
        }

        // 计算下一个频率点
        if (mode == AxisScale::Logarithmic) {
            // 对数模式：freq_next = freq_current * (1 + step_percentage / 100)
            currentFreq += (currentFreq * stepSize / 100.0);
        } else {
            // mode == StepMode::LINEAR
            // 线性模式：freq_next = freq_current + step_value
            currentFreq += stepSize;
        }

        // 如果计算出的下一个点超过了终点，则将下一个点强制设为终点
        if (currentFreq > freqStop) {
            currentFreq = freqStop;
        }
    }
    return pointCount;
}

// 在曲线上插值：xs_Hz 升序（必需）
double CalcUtils::evalAt(double xHz, const QVector<double> &xs_Hz, const QVector<double> &ys,
                         AxisScale xScale, AxisScale yScale, bool *ok)
{
    if (ok)
        *ok = false;
    const int n = xs_Hz.size();
    if (n == 0 || ys.size() != n)
        return 0.0;
    if (xHz <= xs_Hz.front()) {
        if (ok)
            *ok = true;
        return ys.front();
    }
    if (xHz >= xs_Hz.back()) {
        if (ok)
            *ok = true;
        return ys.back();
    }

    // 二分查找
    int l = 0, r = n - 1;
    while (l + 1 < r) {
        int m = (l + r) / 2;
        if (xs_Hz[m] <= xHz)
            l = m;
        else
            r = m;
    }
    double x0 = xs_Hz[l], x1 = xs_Hz[r];
    double y0 = ys[l], y1 = ys[r];

    double X0, X1, X;
    if (!toScaled(x0, xScale, X0) || !toScaled(x1, xScale, X1) || !toScaled(xHz, xScale, X))
        return 0.0;
    double Y0, Y1;
    if (!toScaled(y0, yScale, Y0) || !toScaled(y1, yScale, Y1))
        return 0.0;

    double t = (X1 == X0) ? 0.0 : (X - X0) / (X1 - X0);
    double Y = Y0 + (Y1 - Y0) * t;
    if (ok)
        *ok = true;
    return fromScaled(Y, yScale);
}

QVector<double> CalcUtils::resample(const QVector<double> &xs_Hz, const QVector<double> &ys,
                                    const QVector<double> &gridHz,
                                    AxisScale xScale, AxisScale yScale, bool *ok)
{
    if (ok)
        *ok = false;
    QVector<double> out;
    out.resize(gridHz.size());
    for (int i = 0; i < gridHz.size(); ++i) {
        bool okPoint = false;
        out[i] = evalAt(gridHz[i], xs_Hz, ys, xScale, yScale, &okPoint);
        if (!okPoint) {
            if (ok)
                *ok = false;
        } // 某点失败不影响其他点
    }
    if (ok)
        *ok = true;
    return out;
}

// ====== 单位解析 ======
CalcUtils::UnitSpec CalcUtils::parseUnit(const QString &uIn)
{
    UnitSpec s;
    QString u = uIn.trimmed().toLower();
    s.norm = u;

    // 去掉空格，例如 "dB uV/m"
    u.remove(' ');

    // 功率
    if (u == "w") {
        s.kind = UnitKind::Power;
        s.isDb = false;
        s.scale = 1.0;
        return s;
    }
    if (u == "mw") {
        s.kind = UnitKind::Power;
        s.isDb = false;
        s.scale = 1e-3;
        return s;
    }
    if (u == "uw" || u == "µw") {
        s.kind = UnitKind::Power;
        s.isDb = false;
        s.scale = 1e-6;
        return s;
    }
    if (u == "nw") {
        s.kind = UnitKind::Power;
        s.isDb = false;
        s.scale = 1e-9;
        return s;
    }
    if (u == "dbm") {
        s.kind = UnitKind::Power;
        s.isDb = true;
        return s;
    }
    if (u == "dbw") {
        s.kind = UnitKind::Power;
        s.isDb = true;
        return s;
    }

    // 电压
    if (u == "v") {
        s.kind = UnitKind::Voltage;
        s.isDb = false;
        s.scale = 1.0;
        return s;
    }
    if (u == "mv") {
        s.kind = UnitKind::Voltage;
        s.isDb = false;
        s.scale = 1e-3;
        return s;
    }
    if (u == "uv" || u == "µv") {
        s.kind = UnitKind::Voltage;
        s.isDb = false;
        s.scale = 1e-6;
        return s;
    }
    if (u == "dbv") {
        s.kind = UnitKind::Voltage;
        s.isDb = true;
        return s;
    } // 1 V
    if (u == "dbuv") {
        s.kind = UnitKind::Voltage;
        s.isDb = true;
        return s;
    } // 1 µV
    if (u == "dbmv") {
        s.kind = UnitKind::Voltage;
        s.isDb = true;
        return s;
    } // 1 mV

    // 电流
    if (u == "a") {
        s.kind = UnitKind::Current;
        s.isDb = false;
        s.scale = 1.0;
        return s;
    }
    if (u == "ma") {
        s.kind = UnitKind::Current;
        s.isDb = false;
        s.scale = 1e-3;
        return s;
    }
    if (u == "ua" || u == "µa") {
        s.kind = UnitKind::Current;
        s.isDb = false;
        s.scale = 1e-6;
        return s;
    }
    if (u == "dbua") {
        s.kind = UnitKind::Current;
        s.isDb = true;
        return s;
    } // 1 µA

    // 场强（电场）
    if (u == "v/m" || u == "vpm") {
        s.kind = UnitKind::FieldE;
        s.isDb = false;
        s.scale = 1.0;
        return s;
    }
    if (u == "mv/m" || u == "mvpm") {
        s.kind = UnitKind::FieldE;
        s.isDb = false;
        s.scale = 1e-3;
        return s;
    }
    if (u == "uv/m" || u == "µv/m" || u == "uvpm" || u == "µvpm") {
        s.kind = UnitKind::FieldE;
        s.isDb = false;
        s.scale = 1e-6;
        return s;
    }
    if (u == "dbuv/m" || u == "dbuvpm") {
        s.kind = UnitKind::FieldE;
        s.isDb = true;
        return s;
    } // 1 µV/m

    // 场强（磁场）
    if (u == "a/m" || u == "apm") {
        s.kind = UnitKind::FieldH;
        s.isDb = false;
        s.scale = 1.0;
        return s;
    }
    if (u == "ma/m" || u == "mapm") {
        s.kind = UnitKind::FieldH;
        s.isDb = false;
        s.scale = 1e-3;
        return s;
    }
    if (u == "ua/m" || u == "µa/m" || u == "uapm" || u == "µapm") {
        s.kind = UnitKind::FieldH;
        s.isDb = false;
        s.scale = 1e-6;
        return s;
    }

    // 未识别
    s.kind = UnitKind::Unknown;
    return s;
}

// ====== 驻波比 ======
double CalcUtils::calculateSWR(double fwdW, double revW, bool *ok, QString *errorMessage)
{
    if (ok)
        *ok = false;
    if (errorMessage)
        errorMessage->clear();

    if (!(fwdW > 0.0)) {
        if (errorMessage)
            *errorMessage = "Forward power must be > 0.";
        return std::numeric_limits<double>::quiet_NaN();
    }
    if (revW < 0.0) {
        if (errorMessage)
            *errorMessage = "Reverse power must be >= 0.";
        return std::numeric_limits<double>::quiet_NaN();
    }
    if (revW >= fwdW) {
        if (errorMessage)
            *errorMessage = "Reverse power must be less than forward power.";
        return std::numeric_limits<double>::quiet_NaN();
    }

    // 反射系数 |Γ| = sqrt(Pr / Pf)
    const double gamma = std::sqrt(revW / fwdW);
    // SWR = (1 + |Γ|) / (1 - |Γ|)
    const double swr = (1.0 + gamma) / (1.0 - gamma);

    if (ok)
        *ok = true;
    return swr;
}

static inline double db_to_lin20(double dB)
{
    return std::pow(10.0, dB / 20.0);
}

// 把 dBm 通过 AF 与路径增益/损耗换算成场强 V/m
double CalcUtils::dBm2Vm_general(double dBm, double af_dB, double impedance, double extraPath_dB)
{
    // 等效把路径修正在功率上
    const double dBm_eff = dBm + extraPath_dB;

    // dBm -> V(dBµV)（任意阻抗）
    const double v_dBuV = dBm_eff + 90.0 + 10.0 * std::log10(impedance);

    // 加上天线/探头因子 -> E(dBµV/m)
    const double e_dBuVpm = v_dBuV + af_dB;

    // dBµV/m -> V/m
    return std::pow(10.0, (e_dBuVpm - 120.0) / 20.0);
}

// V/m -> dBm（与 dBm2Vm_general 互逆）
// af_dB: 天线/探头因子，单位 dB/m
// impedance: 系统阻抗（Ω），默认 50
// extraPath_dB: 路径增益/损耗，增益为正、损耗为负
double CalcUtils::vm2dBm_general(double E_Vpm, double af_dB, double impedance, double extraPath_dB)
{
    if (!(impedance > 0.0))
        impedance = 50.0;
    if (!(E_Vpm > 0.0))
        return -std::numeric_limits<double>::infinity(); // 非法场强

    // E(V/m) -> E(dBµV/m)
    const double e_dBuVpm = 20.0 * std::log10(E_Vpm) + 120.0;

    // 去掉 AF：V(dBµV) = E(dBµV/m) - AF(dB/m)
    const double v_dBuV = e_dBuVpm - af_dB;

    // 去掉阻抗常数：dBm_eff = V(dBµV) - [90 + 10*log10(Z)]
    const double dBm_eff = v_dBuV - (90.0 + 10.0 * std::log10(impedance));

    // 源端 dBm：减去路径增益/损耗（增益为正，源端需求更低；损耗为负，源端需求更高）
    return dBm_eff - extraPath_dB;
}

double CalcUtils::dB2Vm(double startVm, double dB)
{
    if (!(startVm > 0.0))
        return 0.0; // 或者返回 NAN
    return startVm * db_to_lin20(dB);
}

double CalcUtils::vm2dB(double firstVm, double secondVm)
{
    if (!(firstVm > 0.0) || !(secondVm > 0.0))
        return std::numeric_limits<double>::quiet_NaN();
    return 20.0 * std::log10(secondVm / firstVm);
}

// 根据参考 dBµA 值和 dB 变化量，计算新的 dBµA 值
double CalcUtils::dB2dBuA(double ref_dBuA, double delta_dB)
{
    if (!(ref_dBuA > -std::numeric_limits<double>::infinity())) {
        if (std::isnan(ref_dBuA) || std::isinf(ref_dBuA)) {
            return std::numeric_limits<double>::quiet_NaN();
        }
    }
    return ref_dBuA + delta_dB;
}

// 计算两个 dBµA 值之间的 dB 差值
double CalcUtils::dBuA2dB(double first_dBuA, double second_dBuA)
{
    // 检查输入有效性 (检查 NaN/Inf)
    if (std::isnan(first_dBuA) || std::isinf(first_dBuA) || std::isnan(second_dBuA) || std::isinf(second_dBuA)) {
        return std::numeric_limits<double>::quiet_NaN();
    }
    return second_dBuA - first_dBuA;
}
